#### ipomopsis ANCOVA example ##### ipo <- read.delim('ecol 563/ipomopsis.txt') out1 <- lm(Fruit~Root+Grazing, data=ipo) library(arm) library(R2jags) x <- ipo$Root y <- ipo$Fruit # create dummy variable for grazing z <- as.numeric(ipo$Grazing)-1 n <- length(y) ipo.data <- list("n", "y", "z", "x") # initial value functions for chains ipo.inits <- function() {list(b0=rnorm(1), b1=rnorm(1), b2=rnorm(1), sigma.y=runif(1) )} # parameters whose posterior distributions should be returned ipo.parms <- c("b0", "b1","b2","sigma.y") # set working director to location of BUGS model setwd("C:/Users/jmweiss/Documents/ecol 563") #### BUGS model ipomodel.txt ##### model { # likelihood for (i in 1:n) { y[i] ~ dnorm(y.hat[i], tau.y) y.hat[i] <- b0 + b1*x[i] + b2*z[i] } # priors b0 ~ dnorm(0, .000001) b1 ~ dnorm(0, .000001) b2 ~ dnorm(0, .000001) # tau.y ~ dgamma(.001,.001) tau.y <- pow(sigma.y,-2) sigma.y ~ dunif(0,10000) } #### end BUGS model ##### # BUGS run ipo.1 <- bugs(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) # JAGS run ipo.1.jags <- jags(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", n.chains=3, n.iter=1000) names(ipo.1) # WinBUGS objects returned by JAGS names(ipo.1.jags$BUGSoutput) # convergence diagnostics ipo.1$summary max(ipo.1$summary[,"Rhat"]) min(ipo.1$summary[,"n.eff"]) # convert to an mcmc object for additional diagnostics bugsfit.mcmc <- as.mcmc.list(ipo.1) jagsfit.mcmc <- as.mcmc(ipo.1.jags) # trace plots of the individual Markov chains xyplot(bugsfit.mcmc) # select only the regression parameters xyplot(bugsfit.mcmc[,c("b0","b1","b2")]) # combined Markov chains in matrix format ipo.1$sims.matrix # combined Markov chains in list format ipo.1$sims.list # individual Markov chains in the order they were obtained ipo.1$sims.array # density plot of posterior distributions for each chain densityplot(bugsfit.mcmc) # percentile credible intervals for parameters colnames(ipo.1$sims.matrix) apply(ipo.1$sims.matrix[,1:4],2,function(x) quantile(x, c(.025,.975))) # compare to frequentist estimates confint(out1) # highest probability density credible intervals # separately by chain HPDinterval(bugsfit.mcmc) # combining the three chains for a single estimate HPDinterval(as.mcmc(ipo.1$sims.matrix)) # compare to percentile intervals apply(ipo.1$sims.matrix[,1:4],2,function(x) quantile(x, c(.025,.975))) # graph of posterior densities for each chain densityplot(bugsfit.mcmc) # graph of posterior densities combining the chains densityplot(as.mcmc(ipo.1$sims.matrix)) # we can also combine chains using rbind and do.call bugsfit.mcmc[[1]] densityplot(as.mcmc(do.call(rbind,bugsfit.mcmc))) #### Bayesian analysis of birds data set ##### birds <- read.csv('birds.csv') names(birds) # removing observations with missing values for response birds.short <- birds[!is.na(birds$S),] dim(birds) dim(birds.short) # simple model of mean richness by year out0.glm <- glm(S~factor(year), data=birds.short, family=poisson) # create variables for WinBUGS y <- birds.short$S # dummy variables for year year2 <- as.numeric(birds.short$year==2006) year3 <- as.numeric(birds.short$year==2007) n <- length(y) # graph of normal priors with decreasing precision curve(dnorm(x,0,sqrt(1/.01)), from=-200, to=200, col=1, ylab='Density') curve(dnorm(x,0,sqrt(1/.001)), add=T, col=2) curve(dnorm(x,0,sqrt(1/.0001)), add=T, col=3) curve(dnorm(x,0,sqrt(1/.00001)), add=T, col=4) #### BUGS model model1.txt ##### model { # likelihood for (i in 1:n) { y[i] ~ dpois(mu.hat[i]) # log link log.mu[i] <- a + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } # priors a ~ dnorm(0, .000001) b1 ~ dnorm(0, .000001) b2 ~ dnorm(0, .000001) } #### end BUGS model ##### # create objects for BUGS run bird.data <- list("y", "year2","year3","n") bird.inits <- function() {list(a=rnorm(1), b1=rnorm(1), b2=rnorm(1))} bird.parms <- c("a", "b1", "b2") # fit models mod0.bugs <- bugs(bird.data, bird.inits, bird.parms, "model1.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=100, debug=T) mod0.bugs <- bugs(bird.data, bird.inits, bird.parms, "model1.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) # check convergence diagnostics # ratio of between to within chain variance max(mod0.bugs$summary[,"Rhat"]) # effective sample size min(mod0.bugs$summary[,"n.eff"]) # examine summaries of posterior distributions mod0.bugs$summary # compare to frequentist estimates coef(out0.glm) # 95% confidence intervals confint(out0.glm) # 95% percentile credible intervals apply(mod0.bugs$sims.matrix,2,function(x) quantile(x, c(.025,.975))) # deviance from Bayesian model: Dbar mod0.bugs$summary["deviance","mean"] # -2 times log-likelihood from frequentist model: Dhat -2*logLik(out0.glm)[1] # number of effective parameters: Dhat - Dbar mod0.bugs$pD # analog of AIC mod0.bugs$DIC #### BUGS model model2.txt: separate intercepts model ##### model { # likelihood for (i in 1:n) { y[i] ~ dpois(mu.hat[i]) # log link log.mu[i] <- a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } # priors for(j in 1:J){ a[j] ~ dnorm(0, .000001) } b1 ~ dnorm(0, .000001) b2 ~ dnorm(0, .000001) } #### end BUGS model ##### # number of patches J <- length(unique(birds.short$patch)) # patch identifier renumbered patch <- as.numeric(factor(birds.short$patch)) # objects for BUGS model bird.data <- list("y", "year2","year3","n", "J", "patch") bird.inits <- function() {list(a=rnorm(J), b1=rnorm(1), b2=rnorm(1))} bird.parms <- c("a", "b1", "b2") # fit model mod1.bugs <- bugs(bird.data, bird.inits, bird.parms, "model2.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=100, debug=T) mod1.bugs <- bugs(bird.data, bird.inits, bird.parms, "model2.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=10000, debug=T) # check convergence diagnostics max(mod1.bugs$summary[,"Rhat"]) min(mod1.bugs$summary[,"n.eff"]) # examine summary of posterior distributions mod1.bugs$summary # compare models with DIC mod1.bugs$DIC mod0.bugs$DIC ##### Random intercepts models ##### #### BUGS model model3.txt ##### model { # likelihood for (i in 1:n) { y[i] ~ dpois(mu.hat[i]) log.mu[i] <- a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } # priors for(j in 1:J){ a[j] ~ dnorm(mu.a, tau.a) } b1 ~ dnorm(0, .000001) b2 ~ dnorm(0, .000001) mu.a ~ dnorm(0, .000001) tau.a <- pow(sigma.a, -2) sigma.a ~ dunif(0, 10000) } #### end BUGS model ##### #### BUGS model model3a.txt ##### model { # likelihood for (i in 1:n) { y[i] ~ dpois(mu.hat[i]) log.mu[i] <- a[patch[i]] + b1*year2[i] + b2*year3[i] mu.hat[i] <- exp(log.mu[i]) } # priors for(j in 1:J){ u0[j] ~ dnorm(0, tau.a) a[j] <- mu.a + u0[j] } b1 ~ dnorm(0, .000001) b2 ~ dnorm(0, .000001) mu.a ~ dnorm(0, .000001) tau.a <- pow(sigma.a, -2) sigma.a ~ dunif(0, 10000) } #### end BUGS model #####